Ejercicio: Datos satelitales en GRASS GIS

Trabajamos con imágenes Sentinel 2 en GRASS GIS

Datos Sentinel 2

Sentinel 2 satellite

• Lanzamiento: Sentinel-2A en 2015, Sentinel-2B en 2017
• Tiempo de revisita: ~5 días
• Cobertura sistemática de áreas terrestres y costeras entre los 84°N y 56°S
• 13 bandas espectrales con resolución espacial de 10 m (VIS y NIR), 20 m (red-edge y SWIR) y 60 m (otras)

Sentinels

ESA - Satélites Copernicus Sentinel. Más información en: https://www.copernicus.eu/en/about-copernicus/infrastructure/discover-our-satellites

Distribución de bandas de Sentinel 2 comparadas con Landsat

Sentinel and Landsat bands

Extensiones para datos Sentinel

• i.sentinel.download: descarga productos Copernicus Sentinel de Copernicus Open Access Hub
• i.sentinel.import: importa datos Sentinel descargados de Copernicus Open Access Hub
• i.sentinel.preproc: importa y realiza corrección atmosférica y topográfica de imágenes S2
• i.sentinel.mask: crea máscaras de nubes y sombras para imágenes S2

Ver Sentinel wiki para más detalles)

Recientemente, se sumaron nuevos miembros en la familia i.sentinel:

• i.sentinel.coverage: comprueba la cobertura de área de las escenas de S1 o S2 seleccionadas

• i.sentinel.parallel.download: descarga imagenes Sentinel en paralelo

• Para conectarse al Copernicus Open Access Hub a través de i.sentinel.download, se necesita ser usuario registrado

• Crear el archivo SENTINEL_SETTING.txt en el directorio $HOME/gisdata/ con el siguiente contenido:

your_username
your_password

Niveles de procesamiento Sentinel 2

• L1C: Reflectancia a tope de atmósfera o Top of Atmosphere (TOA). Disponibles desde el lanzamiento.
• L2A: Reflectancia Superficial o Bottom of Atmosphere (BOA), i.e., los datos han sido corregidos para remover los efectos de la atmósfera. Sólo desde 2019.

Archivo de datos Sentinel

Long Term Archive (LTA)

Todos los productos (1C o 2A) de más de un año son movidos fuera de línea y se requiere un tiempo de espera para ponerlos a disposición del usuario. Esto dificulta la automatización de tareas con productos de más de 12 meses de antigüedad.

Iniciar GRASS GIS, crear nuevo mapset y establecer región computacional

import os

# data directory
homedir = os.path.expanduser('~')

# GRASS GIS database variables
grassdata = os.path.join(homedir, "grassdata")
project = "posgar2007_4_cba"
mapset = "PERMANENT"

import subprocess
import sys

# Ask GRASS GIS where its Python packages are to be able to start it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

# Importar los paquetes python de GRASS
import grass.script as gs
import grass.jupyter as gj

# Iniciar GRASS
session = gj.init(grassdata, location, mapset)

Crear un nuevo mapset sentinel2

# Create a new mapset
g.mapset -c mapset=sentinel2

Definir la región computacionalal radio urbano de Córdoba

# set the computational region to the extent of Cordoba urban area
g.region -p vector=radio_urbano_cba

Búsqueda y descarga de datos S2

Instalar la extensión i.sentinel

# install i.sentinel extension
g.extension extension=i.sentinel

Lista de escenas disponibles que intersectan la región computacional

# explore list of scenes for a certain date range
i.sentinel.download -l \
  settings=$HOME/gisdata/SENTINEL_SETTING.txt \
  start="2020-03-01" \
  end="2020-04-30" \
  producttype=S2MSI2A \
  clouds=30

Lista de escenas disponibles que contienen la región computacional

# filter list of scenes by area_relation=Contains
i.sentinel.download -l \
  settings=$HOME/gisdata/SENTINEL_SETTING.txt \
  start="2020-03-01" \
  end="2020-04-30" \
  producttype=S2MSI2A \
  clouds=30 \
  area_relation=Contains

Descargar la escena seleccionada - NO EJECUTAR

# download the scene that fully contains our region
# i.sentinel.download \
#   settings=$HOME/gisdata/SENTINEL_SETTING.txt \
#   uuid=9a1ea49c-0561-4aa5-ba7a-dc820dc1a316 \
#   output=$HOME/gisdata/s2_data

Como la descarga desde el Copernicus Open Access Hub toma su tiempo, vamos a descargar la escena Sentinel 2 que usaremos y moverla a HOME/gisdata/s2_data

Hagamos una prueba con datos del LTA…

i.sentinel.download -l \
  settings=$HOME/gisdata/SENTINEL_SETTING.txt \
  start="2019-01-01" \
  end="2020-02-28" \
  clouds=30

i.sentinel.download \
  settings=$HOME/gisdata/SENTINEL_SETTING.txt \
  uuid=d4e5df0e-7ead-4407-ba82-d2583be1a6b8 \
  output=$HOME/gisdata/s2_data

Importar datos Sentinel 2 a GRASS GIS

1. Importar con corrección atmosférica: i.sentinel.preproc

Productos nivel 1C

Para obtener un valor de AOD, tenemos 2 opciones:

A. Estimar el valor desde un grafico

B. Descargar un archivo y el valor sera estimado

http://aeronet.gsfc.nasa.gov

Obtener AOD de
http://aeronet.gsfc.nasa.gov

• Estación ARM_Cordoba o Pilar_Cordoba
• Seleccionar fechas de inicio y final
• Seleccionar: Combined file y All points
• Descargar y descomprimir (el archivo final tiene extensión .dubovik)
• Pasar el archivo con la opción aeronet_file

Mapa de elevación

• r.in.srtm.region: importa (y re-proyecta) los mosaicos SRTM que cubren la región computacional, parchea los mosaicos e interpola datos faltantes
• r.in.nasadem: importa (y re-proyecta) los mosaicos de NASADEM que cubren la región computacional y parchea los mosaicos

Si el DEM es más chico que la región computacional, sólo la región cubierta por el DEM será corregida atmosféricamente…

Ejemplo

# enter directory with Sentinel scene and unzip file
cd $HOME/gisdata/s2_data/
unzip $HOME/gisdata/s2_data/name_of_S2_scene

i.sentinel.preproc -atr \
  input_dir=$HOME/gisdata/s2_data/name_of_S2_scene.SAFE \
  elevation=NASADEM \
  aeronet_file=$HOME/gisdata/s2_data/name_of_aeronet_station.dubovik

2. Importar sin corrección atmosférica (as is): i.sentinel.import

Productos nivel 2A

Imprimir información sobre las bandas antes de importarlas

# print bands info before importing
# (1 -proj match, 0 -no proj match)
i.sentinel.import -p input=$HOME/gisdata/s2_data

Importar bandas seleccionadas, recortar y reproyectar al vuelo

# import bands relevant for RGB, NDVI and NDWI
i.sentinel.import -rc \
  input=$HOME/gisdata/s2_data \
  pattern='B(02_1|03_1|04_1|08_1|8A_2|11_2|12_2)0m' \
  extent=region

Listar bandas importadas y revisar metadatos

# list raster maps
g.list type=raster mapset=.

# check metadata of some imported bands
r.info map=T20JLL_20200330T141049_B03_10m
r.info map=T20JLL_20200330T141049_B8A_20m

Balance de colores y composiciones

Asignar grey como paleta de colores

# apply grey color to RGB bands
r.colors \
  map=T20JLL_20200330T141049_B04_10m,T20JLL_20200330T141049_B03_10m,T20JLL_20200330T141049_B02_10m \
  color=grey

Ajuste de colores para una composición RGB color natural

# perform color auto-balancing for RGB bands
i.colors.enhance \
  red=T20JLL_20200330T141049_B04_10m \
  green=T20JLL_20200330T141049_B03_10m \
  blue=T20JLL_20200330T141049_B02_10m \
  strength=95

Mostrar la combinación RGB 432

# display the enhanced RGB combination
d.mon wx0
d.rgb -n \
  red=T20JLL_20200330T141049_B04_10m \
  green=T20JLL_20200330T141049_B03_10m \
  blue=T20JLL_20200330T141049_B02_10m

Tarea

Realizar balance de colores y mostrar combinacion falso color NIR-RED-GREEN

Máscara de nubes y sombras de nubes

Identificar y enmascarar nubes y sus sombras

# identify and mask clouds and clouds shadows: i.sentinel.mask
i.sentinel.mask -s --o \
  blue=T20JLL_20200330T141049_B02_10m \
  green=T20JLL_20200330T141049_B03_10m \
  red=T20JLL_20200330T141049_B04_10m \
  nir=T20JLL_20200330T141049_B08_10m \
  nir8a=T20JLL_20200330T141049_B8A_20m \
  swir11=T20JLL_20200330T141049_B11_20m \
  swir12=T20JLL_20200330T141049_B12_20m \
  cloud_mask=cloud \
  shadow_mask=shadow \
  scale_fac=10000 \
  mtd=$HOME/gisdata/s2_data/S2B_MSIL2A_20200330T141049_N0214_R110_T20JLL_20200330T182252.SAFE/GRANULE/L2A_T20JLL_A016009_20200330T141532/MTD_TL.xml

Visualización de la salida: Nubes y sombras de nubes identificadas con i.sentinel.mask

# display output
d.mon wx0
d.rgb \
  red=T20JLL_20200330T141049_B04_10m \
  green=T20JLL_20200330T141049_B03_10m \
  blue=T20JLL_20200330T141049_B02_10m
d.vect map=cloud fill_color=red
d.vect map=shadow fill_color=blue

Índices de agua y vegetación

Definir región computacional

# set region
g.region -p raster=T20JLL_20200330T141049_B08_10m

Establecer máscara

# set clouds mask
v.patch input=cloud,shadow \
 output=cloud_shadow_mask
r.mask -i vector=cloud_shadow_mask

Estimación de los índices de vegetación

# estimate vegetation indices
i.vi \
  red=T20JLL_20200330T141049_B04_10m \
  nir=T20JLL_20200330T141049_B08_10m \
  output=T20JLL_20200330T141049_NDVI_10m \
  viname=ndvi

Instalar extensión i.wi

# install extension
g.extension extension=i.wi

Estimación de índice de agua

Visualización de los resultados

# estimate water indices and set color palette
i.wi \
  green=T20JLL_20200330T141049_B03_10m \
  nir=T20JLL_20200330T141049_B08_10m \
  output=T20JLL_20200330T141049_NDWI_10m \
  winame=ndwi_mf
r.colors map=T20JLL_20200330T141049_NDWI_10m \
  color=ndwi

Segmentación

Instalar la extensión i.superpixels.slic

# install extension
g.extension extension=i.superpixels.slic

Listar los mapas y crear grupos y subgrupos

# list maps
g.list type=raster pattern="*20200330T141049*" \
  mapset=. output=list.txt

# create groups and subgroups
i.group group=s2 subgroup=s2 file=list.txt

Ejecutar i.superpixels.slic

# run i.superpixels.slic
i.superpixels.slic input=s2 \
  output=superpixels \
  num_pixels=2000

Convertir el resultado a vector

# convert the resulting raster to vector
r.to.vect input=superpixels \
  output=superpixels type=area

Ejecutar i.segment

# run i.segment
i.segment group=s2 output=segments \
  threshold=0.5 minsize=50 memory=500

Convertir el resultado a vector

# convert the resulting raster to vector
r.to.vect input=segments \
  output=segments type=area

Mostrar NDVI junto con las 2 salidas de la segmentación

# display NDVI along with the 2 segmentation outputs
d.mon wx0
d.rast map=T20JLL_20200330T141049_NDVI_10m
d.vect map=superpixels color=yellow fill_color=none
d.vect map=segments color=red fill_color=none

Tarea

Ejecutar cualquiera de los 2 métodos de segmentación con diferentes parámetros y comparar los resultados

Clasificación supervisada

Tarea

• digitalizar áreas de entrenamiento para 3 clases con g.gui.iclass
• guardarlas en un mapa vectorial: training

g.gui.iclass

Clasificación supervisada con Maximum Likelihood

Convertir el vector de áreas de entrenamiento a raster

# convert to raster
v.to.rast input=training output=training \
  use=cat label_column=class

Generar archivos de firma espectral

# obtain signature files
i.gensig trainingmap=training \
  group=s2 subgroup=s2 \
  signaturefile=sig_sentinel

Realizar la clasificación por Maximum Likelihood

# perform ML supervised classification
i.maxlik group=s2 subgroup=s2 \
  signaturefile=sig_sentinel \
  output=sentinel_maxlik

Añadir etiquetas a las clases

# label classes
r.category sentinel_maxlik separator=":" rules=- << EOF
1:vegetation
2:urban
3:bare soil
EOF

Clasificación supervisada con Maximum Likelihood

Clasificación supervisada con Machine Learning

Instalar la extensión r.learn.ml

# install extension
g.extension extension=r.learn.ml

Realizar la clasificación por RF

# perform random forest classification
r.learn.ml trainingmap=training group=s2 \
  output=sentinel_rf n_estimators=300

Añadir etiquetas a las clases

# label classes
r.category sentinel_rf separator=":" rules=- << EOF
1:vegetation
2:urban
3:bare soil
EOF

Clasificación supervisada con Random Forest

Tarea

Comparar los resultados de ambos tipos de clasificación supervisada a través del índice Kappa

Hay un módulo r.kappa

Post-procesamiento y validación

• r.reclass.area para eliminar pequeñas áreas, enmascarar nuevos valores y rellenar los huecos con r.neighbors o r.fillnulls
• convertir la salida en vector y ejecutar v.clean con tool=rmarea
• r.kappa para la validación (idealmente también digitalizar una muestra de prueba)